clc
clear
%%
% This procedure is used to compare the effect of different de-trending
% methods, including multiperiod_EMD,EMD,EEMD,CEEMD,and ICEEMD.
% data0, data1, data2,data3,...,is a multi-period sequence that can be used for
% testing.
% The program has been tested in MatlabR2015a.

% Note:EMD,EEMD,CEEMD,and ICEEMD codes were download from the http://bioingenieria.edu.ar/grupos/ldnlys/metorres/re_inter.htm#Codigos
% If you use this code, refer to the relevant paper.
% This code multiperiod_EMD is written by the author. If used, please refer
% to this paper:LIU Weiqiang, et al. Using a modified empirical mode decomposition 
% method to remove low-frequency trend drift interference in multi-period
% full-waveform induced polarization data.Geophysical journal
% international,In preparation.

% If the code in this article damages your rights and interests, please contact the author in time, 
% we will modify or delete it in the first place.
% Thanks

% LIU Weiqiang
% liuweiqiangedu@126.com , 
% 1228834015@qq.com,
% Chinese Academy of Geological Sciences
% Central South University
% Institute of geology and geophysics, Chinese academy of sciences

%% Download signal and the aimulated trend,data1=data0+trend,data2=data1+[A sudden noise],
load data0.txt                  
load data1.txt
load data2.txt
load data3.txt
load data4.txt
load data5.txt
load data6.txt

load trend0.txt                  
load trend1.txt
load trend2.txt
load trend3.txt
load trend4.txt
load trend5.txt
load trend6.txt

%% Select input data
% SS=data0;
% trend=trend0;
% period=16384;

% SS=data1;
% trend=trend1;
% period=16384;
% 
% SS=data2;
% trend=trend2;
% period=16384;

% SS=data3;
% trend=trend3;
% period=16384;
% 
% SS=data4;
% trend=trend4;
% period=1280*2;
%
% SS=data5;
% trend=trend5;
% period=1280*2;

% SS=data6;
% trend=trend6;
% period=2000*2;


%% An modified empirical modal method for detrending of multiperiod signals 
tic

SS_yuanshi=SS;
for kk=1:3
    [ SS_yuanshi,X_1,X_2] = multiperiod_EMD( SS_yuanshi,period );
end
SS_chuli=SS_yuanshi;
trend_f0=SS-SS_chuli;

toc

figure(1);
subplot(3,1,1);plot([SS,trend_f0],'-');title('The original sequence')
subplot(3,1,2);plot([trend,trend_f0],'-');title('Real trend  and fitting trend')
subplot(3,1,3);plot([SS_chuli],'-');title('The processed sequence')

%%  EMD
tic

IMF1 = emd(SS,'STOP',[0.05,0.5,0.05],'MAXITERATIONS',50);
trend_f1=shiftdim(sum(IMF1(end-1:end,:)));

toc
%%  EEMD
tic

Nstd=0.1;
NR=5;
MaxIter=50;
IMF2 = eemd(SS,Nstd,NR,MaxIter);
trend_f2=shiftdim(sum(IMF2(end-1:end,:)));

toc
%%  CEEMD
tic

Nstd=0.1;
NR=5;
MaxIter=50;
IMF3 = ceemdan(SS,Nstd,NR,MaxIter);
trend_f3=shiftdim(sum(IMF3(end-1:end,:)));

toc
%%  ICEEMD
tic

Nstd=0.1;
NR=5;
MaxIter=50;
SNRFlag=1;
IMF4 = iceemdan(SS,Nstd,NR,MaxIter,SNRFlag);
trend_f4=shiftdim(sum(IMF4(end-1:end,:)));

toc

%%
ts=1/64.*(0:length(SS)-1);
figure(2)
subplot(3,2,1);plot(ts,SS,'k-',ts,trend,'r-');
% ylim([-6,6])
title('Original series')
subplot(3,2,2);plot(ts,trend,'k-',ts,trend_f0,'r.-');
% ylim([-6,6])
title('By multiperiod EMD')
legend('Real trend','Fitting trend')
subplot(3,2,3);plot(ts,trend,'k-',ts,trend_f1,'r.-');
% ylim([-6,6])
title('By EMD')
subplot(3,2,4);plot(ts,trend,'k-',ts,trend_f2,'r.-');
% ylim([-6,6])
title('By EEMD')
subplot(3,2,5);plot(ts,trend,'k-',ts,trend_f3,'r.-');
% ylim([-6,6])
title('By CEEMD')
subplot(3,2,6);plot(ts,trend,'k-',ts,trend_f4,'r.-');
% ylim([-6,6])
title('By ICEEMD')
%%
a=-30;b=30;
a2=-30;b2=30;
figure(5);
subplot(5,4,1);plot(ts,[trend],'k-',ts,trend_f0,'r--');xlim([0,2560]);
% xlabel('Time [s]')
ylabel('\DeltaU [mV]')
title('(a1)')
% ylim([a,b])
subplot(5,4,2);plot(ts,[trend],'k-',ts,trend_f0,'r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(a2)')
% ylim([a,b])
subplot(5,4,3);plot(ts,[trend],'k-',ts,trend_f0,'r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(a3)')
% ylim([a,b])
subplot(5,4,4);plot(ts,[trend],'k-',ts,trend_f0,'r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(a4)')
% ylim([a2,b2])

subplot(5,4,5);plot(ts,[trend],'k-',ts,sum(IMF1(end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
ylabel('\DeltaU [mV]')
title('(b1)')
% ylim([a,b])
subplot(5,4,6);plot(ts,[trend],'k-',ts,sum(IMF1(end-1:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(b2)')
% ylim([a,b])
subplot(5,4,7);plot(ts,[trend],'k-',ts,sum(IMF1(end-2:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(b3)')
% ylim([a,b])
subplot(5,4,8);plot(ts,[trend],'k-',ts,sum(IMF1(end-3:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
% ylim([a2,b2])
title('(b4)')

subplot(5,4,9);plot(ts,[trend],'k-',ts,sum(IMF2(end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
ylabel('\DeltaU [mV]')
title('(c1)')
% ylim([a,b])
subplot(5,4,10);plot(ts,[trend],'k-',ts,sum(IMF2(end-1:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(c2)')
% ylim([a,b])
subplot(5,4,11);plot(ts,[trend],'k-',ts,sum(IMF2(end-2:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(c3)')
% ylim([a,b])
subplot(5,4,12);plot(ts,[trend],'k-',ts,sum(IMF2(end-3:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(c4)')
% ylim([a2,b2])

subplot(5,4,13);plot(ts,[trend],'k-',ts,sum(IMF3(end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
ylabel('\DeltaU [mV]')
title('(d1)')
% ylim([a,b])
subplot(5,4,14);plot(ts,[trend],'k-',ts,sum(IMF3(end-1:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(d2)')
% ylim([a,b])
subplot(5,4,15);plot(ts,[trend],'k-',ts,sum(IMF3(end-2:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(d3)')
% ylim([a,b])
subplot(5,4,16);plot(ts,[trend],'k-',ts,sum(IMF3(end-3:end,:),1)','r--');xlim([0,2560]);
% xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(d4)')
% ylim([a2,b2])

subplot(5,4,17);plot(ts,[trend],'k-',ts,sum(IMF4(end,:),1)','r--');xlim([0,2560]);
xlabel('Time [s]')
ylabel('\DeltaU [mV]')
title('(e1)')
% ylim([a,b])
subplot(5,4,18);plot(ts,[trend],'k-',ts,sum(IMF4(end-1:end,:),1)','r--');xlim([0,2560]);
xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(e2)')
% ylim([a,b])
subplot(5,4,19);plot(ts,[trend],'k-',ts,sum(IMF4(end-2:end,:),1)','r--');xlim([0,2560]);
xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(e3)')
% ylim([a,b])
subplot(5,4,20);plot(ts,[trend],'k-',ts,sum(IMF4(end-3:end,:),1)','r--');xlim([0,2560]);
xlabel('Time [s]')
% ylabel('\DeltaU [mV]')
title('(e4)')
% ylim([a2,b2])


%%
% IMF1=[data5'./1000;SS';IMF1];
% IMF2=[data5'./1000;SS';IMF2];
% IMF3=[SS';IMF3];
% IMF4=[SS';IMF4];

IMF=IMF1;
clf(figure(101))
[aa,bb]=size(IMF);
nn=ceil(aa/3);
for jj=1:aa
    figure(101);
    subplot(nn,3,jj)
    plot(ts,IMF(jj,:),'k-')
%     ylim([-2,2])
    xlim([0,2560])
    ylabel(['IMF' num2str(jj-2)])
end   

IMF=IMF2;
clf(figure(102))
[aa,bb]=size(IMF);
nn=ceil(aa/3);
for jj=1:aa
    figure(102);
    subplot(nn,3,jj)
    plot(ts,IMF(jj,:),'k-')
%     ylim([-2,2])
    xlim([0,2560])
    ylabel(['IMF' num2str(jj-2)])
end    

IMF=IMF3(2:end,:);
clf(figure(103))
[aa,bb]=size(IMF);
nn=ceil(aa/3);
for jj=1:aa
    figure(103);
    subplot(nn,3,jj)
    plot(ts,IMF(jj,:),'k-')
%     ylim([-2,2])
    xlim([0,2560])
    ylabel(['IMF' num2str(jj-2)])
end    

IMF=IMF4(2:end,:);
clf(figure(104))
[aa,bb]=size(IMF);
nn=ceil(aa/3);
for jj=1:aa
    figure(104);
    subplot(nn,3,jj)
    plot(ts,IMF(jj,:),'k-')
%     ylim([-2,2])
    xlim([0,2560])
    ylabel(['IMF' num2str(jj-2)])
end    


%%
% save trend3.txt -ascii trend3
% 
% 
